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Abstract 

A set of Maplev R.3 software routines, for plotting 2D/3D projections of Poincare 
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PROGRAM SUMMARY 



Title of the software package: Poincare. 
Catalogue number: (supplied by Elsevier) 

Software obtainaMe from: CPC Program Library, Queen's University of Belfast, N. Ireland (see applica- 
tion form in this issue) 

Licensing provisions: none 

Operating systems under which the program has been tested: UNIX systems, Macintosh, DOS (AT 386, 
486 and Pentium based) systems, DEC VMS, IBM CMS. 

Programming language used: Maple V Release 3. 

Memory required to execute with typical data: 8 Megabytes. 

No. of lines in distributed program, including On-Line Help, etc.: 1381. 

Keywords: Hamiltonian systems, surface-of-section method, symbolic computing. 

Nature of m,athematical problem 

Computation and plotting of 2D/3D projections of Poincare surfaces-of-section of Hamiltonian systems. 

Methods of solution 

A 4*^^ order Runge-Kutta method with optional stepsize and number of iterations is used. However, it is 
possible to indicate any user-method to be used in the integration scheme. 

Restrictions concerning the complexity of the problem. 

Besides the inherent restrictions of the Runge-Kutta method, this first version of the package does not 
makes use of adaptative stepsize control. 

Typical running time 

It depends strongly on the surface-of-section to be plotted. With a Pentium-90 PC (32 Mb. RAM), fast 
plots usually take from a few seconds to a few minutes. On the other extreme, in an example considered 
in this paper, a surface-of-section with 10,000 points and an energy threshold ~ 10~^ took 35 minutes. 

Unusual features of the program 

This package provides easy-to-use software tools for plottings 2D/3D projections of Poincare surfaces- 
of-section of Hamiltonians systems. The speed at which the returned plots are calculated is adjustable, 
in connection with their accuracy. This feature permits alternatively searching for, say, "first order" 
phenomena at remarkable high speed, or, say, "high order" detailed 2D/3D projections displaying "is- 
lands" and the inner structure of a surface-of-section, as desired. The 2D intersection plane over which 
the surface-of-section is plotted can be any one of the coordinate planes of the phase space, and can be 
shifted in the positive and negative directions. The package also provides routines for setting large sets 
of initial conditions for numerical experiments in seconds. The implementation in a symbolic computing 
environment allows for combined symbolic/numerical studies. 
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LONG WRITE-UP 



Introduction 

The Poincare surface-of-section method has become a popular technique for analyzing 
weakly perturbed Hamiltonian systems 0. People working in correlated areas usually pre- 
pare numerical routines, for instance in Fortran or C, and use the computational environ- 
ment to find the solution curves of the model, and to build the corresponding surfaces-of- 
section plots. In parallel, symbolic systems now present a satisfactory performance when 
computing with hardware floating-point numbers, while offering a comfortable computa- 
tional environment for realizing symbolic mathematical studies and manipulating plots. 

Related to these points, this paper presents a set of software routines, implemented in 
MapleV R.3, for plotting Poincare surfaces-of-section, exploiting the MapleV routines for 
working with hardware floats and with plots. In preparing such a package, our intention 
was: 

• to provide easy-to-use software tools for fast plottings of Poincare sections of Hamil- 
tonian dynamical systems; 

• to make these tools flexible enough to feature sections over any possible 2D/3D 
submanifold of the corresponding phase space (including the time as a possible 
variable) ; 

• to allow for combined symbolic/numerical studies by implementing such software 
tools in a symbolic computing environment. 

The exposition is organized as follows. In Sec.^, the basic ideas concerning the surface- 
of-section method are briefly reviewed. In Sec.|^, a summary of the package's commands 
is presented, followed by a detailed description^ of its most relevant commands, mainly 
poincare, for building the 2D/3D plots of Poincare surface-of-section projections, and 
gin, for generating sets of initial conditions for numerical experiments. Sec.§ contains 
a brief illustration of how the new commands work in three selected examples that can 
also be regarded as simple input/output tests. They are: the Toda lattice the Henon- 
Heiles Hamiltonian 0, and a numerical study, in the context of general relativity, which 
appeared in a recent publication 0. Finally, the Conclusions contain a brief discussion 
about this work and its possible extensions. 

1 The Poincare surface-of-section method 

The Poincare surface-of-section method has become popular in the last decades in con- 
nection with the KAM theorem for weakly perturbed (originally integrable) Hamiltonian 
systems. 

An insight of the ideas underlying this method can be obtained by considering, for 
instance, a conservative system with two degrees of freedom. The physical trajectories 

^Aside from this, the package itself contains an On-Line help in standard Maple format which can be 
viewed as the User's manual for all the routines. 
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lie on the three-dimensional energy surface i?(pi,P2, ^2) = Hq and, if the motion is 
bounded, after a suitable time interval, the solution curves will repeatedly intersect any 
given embedded 2D-plane; for instance, the q2 — constant, plane. If another 

constant of motion involving pi and qi exists (then the system is integrable), it is possible 
to use it to express pi = Pi{qi), i.e., the intersection points will lie on a smooth curve. 
The enclosed area will be an integral invariant and, as time flows, these smooth curves 
will draw surfaces in the phase space. 

For weakly perturbed systems, these surfaces of solution curves of regular motion 
(KAM surfaces) continue to exist for most initial conditions, breaking their topology near 
resonances to form "islands" chains. Within these islands, the topology is broken again 
to other chains and so on. Generally speaking, the KAM surfaces isolate thin layers of 
stochasticity and, as the perturbation strength is increased, transitions between layers 
merge, and the KAM surfaces progressively disappear resulting in complete stochastic 
motion. 

In this manner, the plotting of the intersection points of the solution curves with a 
given 2D-plane (the surface- of-section) permits the study of the existence of non-obvious 
constants of motion (isolating integrals in the context of Hamilton- Jacobi theory), lo- 
cal stability, transition from ordered to stochastic motion, and many other interesting 
properties. We recall that no general procedure for determining the integrability of an 
arbitrary Hamiltonian system, or even the number of such isolating integrals, has yet been 
found. As a consequence, the plotting of these surfaces-of-section plays an important role 
not only in numerical studies, but also in checking the consistency of analytical results 
obtained using perturbative methods. 

In the multidimensional case, the Poincare surface-of-section (PS) has dimensionality 
2N-2. Although the case > 3 presents some subtleties (Arnold diffusion etc.), for 
trajectories which are exactly separable in the (pi, qi) coordinates, the intersection points 
of the solution curves with the corresponding {pi, q^) plane also lie on a smooth curve. We 
will denominate two-dimensional surface- of section (2PS) the plotting of these intersection 
points over a given 2D-plane of the phase space. In the general case the intersection points 
may not lie in such smooth curves, but they still fill an annulus of finite area, and the 
thickness can be related to the nearness concerning exact separability in the related (p^, q^) 
coordinates. 

It is possible to extend the concept of 2D surface-of-section to 3D space-of -section 
(3PS), which will contain 3D projections of the sohition curves. The plotting of the 2PS 
embedded in a related 3PS may give interesting detailed information about the behavior 
of a given model. 

2 The Poincare package 

Basically, the Pomcare package consists of a plotting-command, poinceire, for the 2D/3D 
plotting of the corresponding projections of surfaces-of-section, plus a set of facility-com- 
mands for a quick setup of the Hamilton equations of motion, initial conditions for nu- 
merical experiments, and for the zooming of plots. All plots are built by numerically 
integrating Hamilton's equations for a given set of initial conditions. Once the plot is re- 
alized, all the Maple facilities for handling plots are available. Complementary symbolic 



4 



studies in the same environment can be developed with the commands of the PDEtools 
package m for the analytical solving of systems of ODE's, scalar PDE's (the corresponding 
Hamilton- Jacobi equation), and for changing variables. 

Summary 

A brief review of the commands of the package is as follows^]: 

• hameqs receives a Hamiltonian and returns Hamilton's equations and a list with 
the p's and q's involved; 

• poincare receives a Hamiltonian, a time range and initial conditions, and returns 
a 2D plot of a 2PS; or a 3D plot of 2PS embedded in a related 3PS. In both the 
2D and 3D plots, the 2D plane over which the 2PS is plotted can be any of the 
coordinate planes of the phase space, and can optionally be shifted in the positive 
or negative directions; 

• gin generates a set of lists of initial conditions for plotting a 2PS/3PS, from a given 
incomplete set of fixed values and/or ranges for the p's, g's and the energy. This 
command speeds-up the time-expensive task of setting appropriate initial conditions 
for performing numerical experiments; 

• zoom allows for changing the ranges of the display of a given 2D/3D plot without 
having to recalculate it again, thus saving time and memory resources. 

Description 

A complete description of the Pomcare package's commands is found in the On-Line help. 
Therefore, a detailed description will be given here only for the most relevant routines, 
namely poincare and gin. Some commented input/ output examples can be found in 
Sec.||. 

2.1 Command name: poincare 

Feature: plot 2D/3D projections of the Poincare surface-of-section of a given Hamiltonian 
system. 

Calling sequenc^: 

> poincare(H, t=a. .b, ics, optional_parameters) ; 

Parameters: 

H - any algebraic expression representing the Hamiltonian. 

t=a. .b - t represents the time and a. .b is a numerical range. 

ics - a set of initial conditions for the p's and g's in the form: 

{ [nl ,n2,n3, ...], [...],...}, where [nl ,n2,n3 , . . . ] is a list of numbers 
representing the initial values for [t, pi, p2, . . .pk, ql,q2, . . . ,qk]. 

•^This subsection and the next one may contain some information already presented in the previous 
sections; this was necessary to produce a complete-cut description of the package. 
^In what follows, the input can be recognized by the Maple prompt >. 
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Optional Parameters: 

stepsize=n - indication of a positive number representing the size of the step 

used in the integration scheme. 
iterations=N - indication of a positive integer representing the number of 

iterations used in the integration scheme. 
scene= [xi ,xj] - indication of the variables constituting the 2D-plane of the phase 

space where the 2PS is plotted; it is possible to specify 

ranges for the plot as in: scene=[xi=a. .b,xj=c. .d] 
scene=[xi,xj ,xk] - indication of the 2D-plane and a third variable, xk, to be used as 

cross-variable or third axis in plots of 2PS/3PS's, respectively 
shif t=s - indication of a number representing a positive or negative shift 

of the intersection plane in the plots of 2Ps/3Ps 
method=procedure - indication of a user procedure to be used as integration method. 
3 - to obtain a related 3D plot (3PS). 

Synopsis: 

The poincare command was designed to build either fast but not so accurate, or 
slower, as accurate as desired, 2D/3D projection plots of Poincare sections (see optional 
arguments below). Instead of analytically enforcing the Hamiltonian constraint, the value 
of the energy of each plotted point is checked against the corresponding Hq. As a com- 
plement, another routine, gin (see subsection 3.2), was programmed in order to speed-up 
the preparation of initial conditions for the numerical experiments. 

The returned plots can be manipulated using all the standard Maple facilities (icon 
tool-bar) for 2D/3D plots such as reorientation, perspective, etc., and using the zoom 
command, also part of this package, for zooming in/out the interesting regions. These 
facilities usually permit a detailed visual distinction between the KAM surfaces and the 
layers of stochasticity for, say, typical weakly perturbed systems. 

In addition to the returned plot, some relevant information related to each solution 
curve is displayed on the screen during the calculations. This information is: 

• the initial value of the Hamiltonian, p's and q's for the curve; 

• the number of intersection points found in the given time interval (2D plots); 

• the maximum percentile "energy-deviation" of the intersection points. 

It is useful to know the number of intersection points since it may be an indication of how 
appropriate is the indicated time interval. Concerning the percentile energy-deviation]^, 
it is calculated as ^o-H^pomt) -j^gg note that all numerical algorithms will lead to 

values of H different from the initial Hq, specially in the case of optional fast plottings. 
Percentile deviations below 10^*^ are displayed as 0. Also, the deviation is calculated 
for each intersection point, but only the greatest value is displayed. This information 
will give an idea of the accuracy of the plot, and the user will have the option of either 
reentering the instruction looking for a more accurate/slower plot or using his/her own 
numerical integration algorithm (see below). In typical situations, the smooth curves can 
be recognized even with energy-deviations of the order of Hq/1Q. 

^When Hq — 0, just the maximum absolute deviation is displayed. 
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The arguments 

The first argument of a poincare calling is the Hamiltonian. Some handy conventions 
were adopted to represent the g's and p's. These are: all p's and g's must appear as pn 
or qn where n is a positive integer, as in pi, p2, etc., and the time dependence need not 
be explicit, as in pn or qn instead of qn{t) or pn{t). The Hamiltonian is assumed to be 
time independent and the number of degrees of freedom^ is expected to be < 10. 

The solution curves are calculated within a given time interval specified in the second 
argument as t=a. .b. When this time range leads to no intersection points, one can use 
the 3 option, see how far the intersection plane is from the trajectories, and reenter the 
instruction shifting the intersection plane using the shift option (see below). 

The third argument, a set, may have any number of lists of initial conditions for 
the time, p's and q's, corresponding or not to the same initial value of H; they can be 
generated by the user or by the gin command, as explained in subsection p.2| . The initial 
conditions must be inside a set structure as in {[nl,n2,n3, . . .] , [...],...}, where 
[nl,n2,n3, . . .] is a list of numbers representing the initial values for [t, pi, . . .pk, 
ql , . . . qk] . 

The optional arguments 

The optional arguments can be given alone or in conjunction and in any order. 

By default, the step interval is (b-a)/20, where a. .b is the range for t. This can be 
changed by giving the extra argument stepsize = n, where n is a positive number. As 
the stepsize is decreased, the accuracy and the smoothness of the integral curves (as well 
as the time consumed in the calculations) will increase. 

The default numerical algorithm used in the integration scheme is basically the 4*^^ 
order Runge-Kutta of the DEtools package, but this can be changed by giving the extra 
argument method = usermethod, where usermethod should be a numerical integration 
algorithm (see the help pages of the DEtools package). Also, when requesting a plot, 
the numerical algorithm can be iterated, as many times as desired, by giving the extra 
argument iterations = N (default: iterations = 1). 

The default scene for the plots is: the {pi^qi) plane, at g2 = for the 2PS, or the 
{pi,qi,q2) 3D-submanifold, for the plot of a 2PS embedded in a 3PS, when the 3 option 
is indicated. The intersection points constituting the 2PS are obtained by looking for 
the sign change of a "third" coordinate, here denominated the cross-variable, by default 
q2- The default ranges for the display of a plot are such as to include all the calculated 
intersection points in the case of 2D plots; or all calculated pieces of projections of solution 
curves plus the intersection 2D-plane, in the case of 3D plots. 

All these defaults for the scene can also be changed. First of all, the intersection plane 
over which the 2PS is plotted can be shifted in the positive or negative directions by 
indicating shift=s (a real number) as an additional argument. Concerning the 2D/3D 
submanifolds or the cross-variable, they can be changed by giving the extra argument 
scene= [xl ,x2] or scene= [xl ,x2,x3] , with or without extra ranges (for displaying a 
plot of a particular region), as in scene=[xl=a. .b, . . .]. The 2PS will then be plotted 
over the plane formed by the first two variables appearing in the right-hand-side, and 
X3, when given, will be the cross-variable; or the third axis when the 3 option is given as 

^To avoid useless large code this number was restricted to 10 (20-dimensional phase space), but this 
restriction can easily be removed. 
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argument tooQ. When ranges are indicated, it is still possible to zoom in/out the resulting 
plot, up to the default ranges mentioned above, by using the zoom command. 

2.2 Command name: gin 

Feature: generates a set of lists of initial conditions satisfying the Hamiltonian constraint. 
Calling sequence: 

> gin(H,ic,N); 

Parameters: 

H - any algebraic expression representing the Hamiltonian. 

ic - a set of single initial conditions for the time t, p's, g's and the energy, in the form: 
{t=Xl,pl=X2, . . ,qn=Xk,energy=X2n+2)}, where the Xi are numbers or ranges of 
numbers, representing initial values for the time, p's, g's, and the energy. 

N - a positive integer representing the number of desired lists of initial conditions 

Synopsis: 

The numerical study of a given Hamiltonian system usually requires a lot of suitable 
lists of initial conditions (IC's) satisfying the Hamiltonian constraint (see examples in 
SecH). The gin command was thus proposed as a tool for the fast generation of a set of 
such lists, with the appropriate syntax as required by the poincare command explained 
above. 

The first argument received by gin is the Hamiltonian. As second argument, a single 
specification of IC's for all but one of (t,p's,q's, Hq) must be given. If IC's are specified 
for all these variables, the command just checks the values against H for consistency. The 
IC for a given variable may be a fixed number or a numerical range. 

The third argument is the number of requested lists of IC's. Face the request of 
say IC's giving ranges or fixed numbers for all variables but one, say x, the routine 
proceeds as follows. The received IC's are separated into numerical ranges and fixed 
numbers. Each range is then uniformly divided into N numbers, and lists are built by 
taking one number, from each divided range, together with the received fixed numbers 
for the other variables. These lists are sequentially introduced in H, resulting in A^ 
algebraic equations for x. Each equation is then numerically solved, and a set of A^ lists 
of complete IC's satisfying the Hamiltonian constraint, with values inside the given ranges, 
is returned. Remarkably, though gin is a very simple command, it plays a fundamental 
role in speeding-up the numerical studies. 

3 Examples 

This section contains a brief illustration of how the routines here presented work in some 
selected examples: the Toda lattice, the Henon-Heiles Hamiltonian, and a numerical study, 

^It is possible to indicate the time "t" as the third variable, in which case a convenient mouse- 
manipulation of the 3D-plot can display the projection of the curves over each (g^, qj) plane. This may 
be useful to study the bounded/unbounded properties of a given potential. Furthermore, in the case of a 
system with 3 degrees of freedom, the use of the "3" option with sccne=[ql,q2,q3] will render the 3D-plot 
of the physical trajectory. 
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in the context of general relativity, which appeared in a recent publication 0]. Although 
the idea is just to test the obtainment of some of the related PSs using the package's 
commands, we also included, for completeness, brief introductions to each example as 
well as short comments about the resulting PS plots. 

3.1 The Toda lattice 

Our first example is the three-particle Toda lattice[Q]. This model describes the motion of 
three particles of equal mass moving on a ring with exponentially decreasing interactions. 
Early numerical experiments]^], followed by analytical results^, showed that the related 
dynamical system is integrable. After some manipulations, the original Hamiltonian can 
be written as: 

2 24 \ / 8 

A fast plot of the surfaces-of-section, projected over the {p, q) planes, with just one initial 
condition, can be obtained via: 

> H, t=-150..150, {[0, .1,1.4, .1,0]}; # ics: t=0,pl= . 1 ,p2=l . 4,ql= . 1 ,q2=0 

> poincareC, stepsize= . 05 , iterations=5) ; 

> poincare ( " " , stepsize= . 05 , iterations=5 , scene= [p2 , q2] ) ; 



Fig.l.a. 2PS over the 52=0 plane, with 127 inter- Fig.l.b. 2PS over the qi=0 plane with 146 inter- 
section points lying on smooth curves. -ffo=0-99, section points. Time: 76s. 
iJ-deviation=5xlO-6%. Time: 73s. 

Generally speaking, it is interesting to have the 2PS projected over all the {pi, qi) planes 
of the phase space, since the existence of smooth patterns in all of them is an indication 
of the possible integrability of the system. 

A Poincare space-of-section corresponding to Fig.l.a can be manipulated using the 
mouse to obtain the following illustrative perspective^: 

®In what follows, the values of 6 and </> mentioned in the figures correspond to the Maple values for 
the plots. 
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> poincare(H,t=-100. .100, {[0, .1,1.4, . I,0]},stepsize=. l,iterations=4, 

> scene=[pl=-1.5. . 1 .5,ql=-l .5 . . 1 . 5, q2=-l .2 . .1.3] ,3) ; 



Fig. 2. a. 3PS projected at 6=—20, (f>=75, showing Fig.2.b. A plane projection of the 3PS shows how 
a KAM surface of regular trajectories. Time: 41s. the intersection points are joined outside the 2PS. 
Another indication of the integrabihty of the system is that regular curves exist what- 
ever the value of Hq. As an example of this, a surface-of-section (one solution curve) and 
a related 3PS, at Hq = 256, can be built as follows: 

> ics2 := gin(H,{t=0,p2=22,ql=0,q2=0,energy=256},l) : 

> poincare(H,t=-50 . .50,ics2,stepsize=.005,iterations=4,scene=[p2,q2] ) ; 

> poincare(H,t=0 . .20,ics2,stepsize=.01,iterations=4,scene=[p2,q2,ql] ,3) ; 



Fig. 3. a Smooth curves on the 2PS, ^1=0 plane. Fig.3.b. 3PS corresponding to Fig. 3. a., 6'=100, 
i?o=256, i?-deviation=3xl0~^ %, 342 points. 0=40, displaying a KAM surface constituted by 
Time: 265s. just one regular curve. Time: 40s. 

3.2 Henon-Heiles Hamiltonian 

The Henon-Heiles Hamiltonian produces one of the most famous and studied surfaces-of- 
section. The corresponding Hamiltonian can be obtained by expanding the Toda Hamil- 
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tonian to cubic terms in qi and q2, and is given by: 




(2) 



This model can be related to the motion of a star in a cylindrically symmetric gravitational 
galactic potential, is not integrable, and the phase space is bounded only if the energy 
is less than 1/6. Well-studied surfaces-of-section, presented in several treatises of chaos, 
with Hq equal to 1/24, 1/18, 1/12, 1/8, 1/7, and 1/6, are obtained here using the gin 
and poincare commands as follows. To start with, six sets, related to each value of Hq 
respectively, with three different initial conditions each, are generated via: 



> for h in [1/24,1/18,1/12,1/8,1/7,1/6] do 

> ics[h] := gin(H,{t=0,p2=0. l,q2=-0.2. .0.2,ql=-0.2. .-0. l,energy=h},3) 



After that, surfaces-of-section with around 550 points, calculated in approximately 5 
minutes each, with percentile if-deviations ~ 10~^ %, can be obtained via: 

> for h in [1/24,1/18,1/12,1/8,1/7,1/6] do 

> poincare (H,t=-300. .300, ics[h] , 



> od: 



> 



stepsize=.05,iterations=3,scene=[p2=-.5. .0.5,q2=-.5. .0.5] ) ; 



> od: 



Fig.4.a. Ho = 1/24 



Fig.4.b. Ho = 1/18 



Fig.4.c. Ho = 1/12 



Fig.4.d. Ho = 1/8 



Fig.4.e. Ho = 1/7 



Fig.4.f. Ho = 1/6 
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The figures above reflect the progressive disintegration of the KAM surfaces, occurring 
with the increase of Hq up to 1/6. In the plots for Hq equal to 1/24 and 1/12, invariant 
curves apparently exist everywhere, but this is not strictly correct. In fact, the model 
is not integrable, as is reflected by the sequence of figures, and thin resonance layers of 
stochasticity are densely distributed throughout the 2PS, even for small i^o- 

Another interesting parameter is given by the relevant number of calculations involved 
in the building of each of Figs.4.(a,b,c,d,e,f): 600/0.05=12,000 points, iterated 4x3=12 
times each; that is, 144,000 calculations in five minutes, or 480 calculations per second]^. 

Concluding this example, an instructive test is to compare the regularity of the physi- 
cal trajectories of Fig.4.a {Hq = 1/24) with that of the trajectories of Fig.4.f. {Hq = 1/6). 
These six displays are obtained by plotting "3PS's" over {qi,q2,'t) and changing the per- 
spective, appropriately, using the mouse: 

> for ic in [ics [1/24] [1 . . 3] , ics [1/6] [1 . . 3] ] do 

> poincare(H, t=-100 . . 100 , {ic} , scene= [ql ,q2,t] , steps ize= . 05 , iterations=3,3) : 

> od: 



Fig.5.a. i/o = l/24, ics[l/24][l], Fig.5.b. i7o=l/24, ics[l/24][2], Fig.5.c. Ho = l/2i, ics[l/24][3], 
H-dev.^2AQ-'' %, time: 124s. i7-dev.=2xl0"'^ %, time: 149s. i7-dcv.=3xl0"'^ %, time: 159s. 



Fig.S.d. ffo=l/6, ics[l/6][l], H- Fig.S.e. Hn=l/6, ics[l/6][2], H- Fig.5.e. Ho=l/6, ics[l/6][2], H- 
dev.=4xl0"'^ %, time: 169s. dev.=5xl0"'^ %, time: 178s. dev.=6xl0"'^ %, time: 191s. 

^The number of digits carried in floats is 16; also, the true calculations per second is greater than 480: 
part of the flve minutes were dedicated to determine the intersection points and the maximum percentile 
deviation. 
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The contrast between Figs.5.(a,b,c) and Figs.5.(d,e,f) emphasizes the connection be- 
tween the smoothness of the curves over a surface-of-section and the regularity of the 
corresponding physical trajectories. 

3.3 An example from general relativity 

Chaos in General Relativity is an issue of rich debate. Since the pioneering work of Beli- 
inski, Khalatnikov and Lifshitz[0], who discussed chaotic behavior in anisotropic Bianchi 
IX models, much advance has been made. Francisco and Matsas[|l, studying that model, 
originally noted that the Liapunov exponents tend to zero in the numerical experiments, 
due to the choice made for the time variable. This result has been confirmed by several 
authors^. In this context, the Poincare surfaces-of-section play a crucial role, since the 
destruction of KAM surfaces does not depend on the time variable and constitutes the 
signature of chaos. 

Bearing the above remarks in mind, we examine here a numerical study appearing in 
a recent work by Calzetta and El Hasi0]. The authors developed a perturbative study 
of the influence of the scalar radiation field on the expansion of the universe in the early 
stages of inflation. They performed numerical experiments to exhibit chaotic behavior 
indicated by the destruction of tori structures, formation of cantori, and Arnold diffusion. 
The Hamiltonian for the model is given by: 

H=^ {-pI -ql + 2Aqt + pl + ql + w?qlql) = (3) 

Such an expression describes closed homogeneous and isotropic universes with a cosmo- 
logical constant. A, playing the role of the infiaton. The degrees of freedom of the model 
are the "radius of the universe", gi, and the conformally scaled radiation field, 
represents its mass, and pi, p2 are the conjugated momenta. 

Taking convenient values m = 0.65 and A = 0.125 and choosing the q2 = 0, {qi,pi) 
plane as in 0, it is possible to reproduce the relevant 2PS there displayed as Fig.2 as 
follows. First, an appropriate set of one hundred lists of initial conditions, satisfying 
Ho = and in accordance with the initial values for pi, qi indicated in that paper, is 
generated via: 

> icl := gin(H,{t=0,pl=-.2. .-0.7,ql=0,q2=0,energy=0},15) : 

> ic2 := gin(H,{t=0,pl=-.7. .-0.812,ql=0,q2=0,energy=0},85) : 

> ics:= icl union ic2: 

> poincare (H , t=0 . . 300 , ics , 

> stepsize= . 1 , iterations=3, scene= [ql=-l . 5 . . 1 . 5 ,pl=-l . . 1 ,q2=-l ..!]); 
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Fig. 6 The presence of smooth curves related to KAM surfaces, as weh as a region of broken 
tori. More than 10,000 intersection points, absolute if-dev.«10~*. Time: 35 minutes. 

Another good test is given by the plotting of a 3PS over (pi, qi, g'2), for initial conditions 
very close to the critical point Pc, where pi — P2 — Qi — Q2 — 0: 

> ics := gin(H,{t=0,pl=0. .-0.31,ql=0,q2=0,energy=0},40) ; 

> poincare (H , t=0 . . 20 , ics , 

> stepsize=. l,scene=[ql=-.3. .0.3,pl=-.3. .0.3,q2=-.l. .0.1] ,3); 



Fig.T.a. 6=— 160, 0=50. Regular circles close to 
Pc, in agreement with complex eigenvalues for the 
linearized system. Time: 61s. 



Fig.7.b. 6=—20, (/>=18. The action of higher order 
terms in splitting- up the KAM surfaces. Absolute 
//■-deviation < 10~^. 



14 



4 Conclusions 



This work presented a set of software-tools for the plotting of Poincare surfaces-of-section 
of Hamiltonian systems. This method is a valuable tool for studying dynamical systems, 
since it conveys relevant information on the dynamics, in a simple manner. Due to the 
characteristic flexibility of symbolic programming languages, the main result was an easy- 
to-use package of commands permitting reasonably fast and varied numerical studies of 
Hamiltonian systems, in a general purpose symbolic computing environment. 

An important remark, taking into account that this software was written for realizing 
intensive numerical studies, is that great emphasis was put in the interactive character of 
the package. That is, the user is given the possibility of alternatively searching for "first 
order" phenomena at remarkably high speed, obtaining draft Poincare sections in just a 
few seconds; or "high order" detailed 2D/3D projections, displaying "islands" and the 
inner structure of a PS, as desired. 

On the other hand, this is a first version of the package and, as such, it docs not 
make use of the theory of Liapunov coefficients, does not discuss the determination of the 
analytical Poincare mappings, does not feature options related to the structure of Arnold's 
web, and is not designed to study the interesting case of dissipative Hamiltonians and 
correlated phenomena (strange attractors etc.). All these topics are possible extensions 
of this work and we expect to report related work in the near future. 
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